Analysis of the Disorder-Induced Zero Bias Anomaly in the Anderson-Hubbard Model 
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Using a combination of numerical and analytical calculations, we study the disorder-induced 
zero bias anomaly (ZBA) in the density of states of strongly-correlated systems modeled by the 
two dimensional Anderson-Hubbard model. We find that the ZBA comes from the response of the 
nonlocal inelastic self-energy to the disorder potential, a result which has implications for theoretical 
approaches that retain only the local self-energy. Using an approximate analytic form for the self- 
energy, we derive an expression for the density of states of the two-site Anderson-Hubbard model. 
Our formalism reproduces the essential features of the ZBA, namely that the width is proportional 
to the hopping amplitude t and is independent of the interaction strength and disorder potential. 



I. INTRODUCTION. 

The Anderson-Hubbard model (AHM) is the simplest 
model that describes strongly-correlated electrons in a 
disordered lattice. The AHM is widely used, for exam- 
ple, to describe doped transition metal oxides, where the 
electronic properties are affected by both a strong lo- 
cal Coulomb repulsion and doping-related disorder .i The 
AHM is also relevant to cold atomic gases in random op- 
tical latticesr'"— and there has been recent interest in the 
AHM as a model interacting system that exhibits An- 
derson localizationj^"— The physics of the AHM is de- 
termined by dimensionality, by filling, and by the three 
energy scales: the kinetic energy t, the on-site Coulomb 
repulsion C7, and the disorder strength A. When A = 0, 
this model reduces to the well-known Hubbard model 
which, despite its simplicity, has only been solved exactly 
in the limits of one^^ and infinite dimensions'^. 

In the Hubbard model, the interesting physics arises 
from a competition between t, which tends to delocalize 
electrons, and U , which tends to localize electrons. When 
the lattice is half filled (i.e. when there is one electron per 
site) , a sufficiently large U can generate a Mott insulating 
phase. The Mott transition occurs at a critical Uc that 
depends on the details of the lattice. Much of the Hub- 
bard model research in the past few decades has revolved 
around strong correlation effects slightly away from the 
Mott insulating phase, which is achieved either by taking 
U less than Uc or by doping away from half filling. One 
of the important ideas to come out of the Hubbard model 
is that the low energy physics of the strongly correlated 
metal phase near the Mott transition is governed by an 
effective interaction J ^ t^/U^ 

Contrary to this, recent exact diagonalization and 
quantum Monte Carlo studies of the two-dimensional 
Anderson-Hubbard model have found that a zero bias 
anomaly (ZBA) of width t forms in the density of states 
(DOS)-^^ The ZBA appears as a V-shaped dip in the 
DOS at the Fermi energy e^, as shown in Fig. [T] While 
it is not surprising that disorder might introduce another 
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FIG. 1: (color online) Density of states for electron densities 
n = 1 (half filling) and n — 0.8, showing the V-shaped zero 
bias anomaly ai ef- Results are for exact diagonalization of 
12 site lattices, and are averaged over 1000 disorder configura- 
tions. Model parameters are A = 20t and U = 8t throughout 
this work, unless stated otherwise. 



low energy scale other than /U , it is surprising that this 
new scale is independent of both U and A. 

It is worth emphasizing that the observed ZBA is not 
explained by the conventional Altshuler-Aronov theory 
of weakly correlated metals. In Altshuler-Aronov the- 
ory, the magnitude of the ZBA depends inversely on a 
dimension-dependent power of the Fermi velocityf22, while 
the AHM ZBA grows linearly with the Fermi velocity 
(which is approximately 2t). 

The physics of this ZBA is subtle, and is not 
captured by most approximations. The Hartree- 
Fock approximation^i^^"— yields a V-shaped zero bias 
anomaly when magnetic moments are allowed to form^^ 
and has a low-energy soft gap that is apparently associ- 
ated with a multi-valley energy landscape. However, 
the width of the ZBA grows with U , suggesting that 
the physics of the ZBA is different than that found by 
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exact diagonalization. Furthermore, the evidence for a 
soft gap in exact diagonahzation calculations is less well 
established,— and it is possible that quantum fluctua- 
tions fill in the soft gap. Another common approxima- 
tion, dynamical mean field theory (DMFT),i^ii^ii^i22,-^ 
includes strong correlation physics, but has not found a 
ZBA at all. It has been argued'^'* that this is because 
of nonlocal contributions to the self-energy neglected in 
these calculations. Recent analytical studies of the two 
site AHM do find a ZBA with qualitative features that 
are consistent with exact diagonalization. These calcu- 
lations interpret the ZBA in terms of level repulsion be- 
tween many-body eigenstatesj^^""— and demonstrate how 
strong correlations can generate a kinetic energy driven 
ZBA. While these studies are instructive, it is difficult to 
connect them to the more usual language of many-body 
self-energies in interacting systems. 

In this article, we show how the ZBA arises from the 
response of the inelastic self-energy to the disorder po- 
tential, using an approach that is loosely based on one 
used by Abrahams et ali^ to study the ZBA in weakly- 
correlated metals. We restrict ourselves to two dimen- 
sions, where the existence of the ZBA is well established, 
and work in the limit of strong disorder. In Sec. Ill A| we 
show that the ZBA comes from nonlocal contributions to 
the local density of states, establishing (i) that the ZBA 
is not a remnant of the Mott gap and (ii) that approxima- 
tions such as Hartree-Fock and DMFT (which retain only 
the local self-energy) are missing key nonlocal physics. In 
Sec. Ill B1 we discuss an approximate self-energy, based on 
equation-of-motion calculations?^ which highlights the 
role of nonlocal spin and charge correlations. We show 
numerically that this approximation works well for large 
disorder, and then derive in Sec. Ill Cl an approximate ex- 
pression for the density of states (DOS) based on this 
self-energy. We find that the energy t appears as the 
natural energy scale for the ZBA. The results are sum- 
marized in Sec. IIIII 



II. CALCULATIONS 

Before we proceed with the calculations, we emphasize 
a significant difference between weakly and strongly cor- 
related systems that affects our analysis. In the atomic 
limit, obtained by setting t = 0, the DOS is a sum of the 
local spectrum at each atomic site. For noninteracting 
systems, each local spectrum has a single resonance at 
the orbital energy e of that site. However, for strongly 
correlated systems, there are two resonances, at e and 
e + U, which we term the lower Hubbard orbital (LHO) 
and upper Hubbard orbital (UHO) respectively. These 
energies correspond to transitions in which an electron is 
added to a site that is initially empty (LHO) or singly 
occupied (UHO). The LHO and UHO are precursors of 
the lower and upper Hubbard bands that form when t is 
nonzero. 

The calculations in this work are based on an expan- 



sion around the atomic limit and are appropriate for 
the strong disorder case. By strong disorder, we mean 
A/2z ^ t, where z is the coordination number of the 
lattice and A/2z is of the order of the average level spac- 
ing of the z sites adjacent to any site in the lattice, and 
the factor of 2 is because there is an LHO and a UHO at 
each site. In this limit, the local spectrum at a particular 
lattice site is dominated by 2(z-|-l) resonances associated 
with the site and its z nearest neighborSf^^ 

A. Analysis of Numerical Results 

In this section, we develop a framework that explicitly 
shows the role of local and nonlocal correlations in the 
DOS. We then use this framework to analyze the results 
of numerical exact diagonalization calculations for the 
AHM. We begin with a brief description of the exact 
diagonalization calculations. 

The AHM Hamiltonian is 

^ = ^ U]c\^Cj„ + ^ (eiUj -I- U fii^fiii) , (1) 

where Uj = —t for nearest- neighbor sites i and j, and 
is zero otherwise; and the annihilation and 

number operators for lattice site i and spin a, and is 
the energy of the orbital at site i. Disorder is introduced 
by choosing from a uniform distribution Ci G tI' 
The AHM can be solved exactly for small clusters. 
For our numerical work, we use a standard Lanczos 
method^i to find the ground states of two-dimensional 
iV-site {N ~ 10, 12) clusters with periodic boundary 
conditions, and then use a block-recursion method to find 
the full nonlocal Green's function Gij (w) for the lattice.— 
The DOS is 

p{E)^-^ln^{Y,Gu{E)), (2) 

i 

where (...) indicates an average over disorder config- 
urations at fixed chemical potential. Examples of the 
disorder-averaged DOS are shown in Fig. [T] 

The goal of this section is to relate the DOS to two 
physically interesting quantities, the local inelastic self- 
energy Tiii{uj) and the nonlocal hybridization function 
Ai(cj). For a given disorder configuration, the inelastic 
self energy is 

S(c^) -Go(c^)-i -G(^)-\ (3) 

where (. . .)^^ is a matrix inverse, Go(w) is the noninter- 
acting Green's function for the same disorder configura- 
tion as G(lli), and Tni{u) is a diagonal matrix element 
of S(a;) (bold symbols indicate matrices in the space of 
lattice sites). The hybridization function is then defined 

by 

Gu{uj) - [c^ - e, - - A,(c^)]-i, (4) 
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where Gii{u) is the local Green's function at site i, and 
T,n is a diagonal matrix element of Both 
and Ai(w) can be extracted from our numerical calcu- 
lations: Eq. ([3]) gives Eii(a;), and then Eq. (jH) can be 
inverted to find Ai{uj). 

In the following analysis, we derive a formal expression 
for p{E) in terms of and Ai{uj). Our starting point 

is Eq. ([2]), with Gii{uj) given by Eq. It is clear from 
these two equations that p{E) depends directly on T,ii{uj) 
and Ai(a;), and the main issue we face in our derivation 
is how to perform the disorder average in Eq. We do 
this in two steps: first, we take a partial disorder average 
of Sii(w) and Ai(w) over ej for j ^ i and for fixed e^; 
second, we average Ge. [uj) over e^. As a result of the first 
averaging process. 



(5) 



where Si{uj) = Tju{uj) + Ai(a;). This gives the average 
self-energy of all sites with energy e. Then 



(6) 



Equation ([5]) is the main approximation made in our 
derivation, and we check below that we do not lose the 
physics of the ZBA as a result of it. The next step is to 
average Ge(w) over the local site energy. 

To perform this average, we expand S'e(a;) about an 
energy E near e^?, by analogy to what is done in Fermi 
liquid theory. In making this expansion, we consider two 
categories of site: (i) sites with e ^ E (LHO near E) 
and (ii) sites with e + U ^ E (UHO near E) . Sites with 
neither e nor e + U near E do not contribute to the DOS 
at E and are not included in our calculations. For cases 
(i) and (n) 

5,H « S^{E) + {e-E)d,S,{E)^^^ 

+ {lo - E)d^S^{Lu)^=E (7) 

where i? = £' for case (i) and E = E — U for case (ii). 
Then the local Green's function for site energy e is 

Z 



E-{e- E)/m* - Z[S^{E) ~ U] 



(8) 



with Z = [l-d^_S-^{u;)]zlE: = Z[l + d,S,iE)] 

and U ~ E — E. The final term in the denominator, 
S^{E) — U, vanishes identically in the atomic limit (Ap- 
pendixl^l. Near the atomic limit, S'e(w) is complex, with 
small real and imaginary parts that shift and broaden the 
orbital energies. We show in Sec. Ill Bl that the imaginary 
part of Sc{E), which results from disorder averaging, is 
of order zt'^/A. 

Because the imaginary part of Se{E) is small, the av- 
erage of Ge {E) over e is easily done (Appendix [BJ , and 
we obtain the DOS 
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FIG. 2: (color online) Origin of the zero bias anomaly. Data 
are for n = 1, U = 8, A = 2Qt, t = 1, unless otherwise 
indicated, (a) Ee^ [E) for E = ef = U /2 versus site en- 
ergy ei. Data is shown for sites with ei + U ~ E (left) and 
ei ~ E (right). Lines are quadratic least-squares fits to the 
data in each region and are used to determine (?eEe(i5); (b) 
9eEe(_B)|g^-g (squares) and 9eAe(iJ)|^^-g (circles) ior E = E 
(solid symbols) and E = E — U (empty symbols), (c) The 
approximate DOS, calculated using the results from (b) in 
Eq. ((9} (solid circles), the exact DOS (solid line), and the 
DOS for the approximate self-energy (I14p (empty circles) . (d) 
Results of a similar analysis for A = 8. 



where we have adopted the convenient notation 



55lho(-B) 



Red,S,iE)U^E 
Red,S,iE)U^E-u- 



The two terms in the sum in Eq. (|9]) give the partial 
DOS for the LHO and UHO. For each term, there are 
two distinct contributions: the first, dfT,^(E), includes 
local Mott physics; the second, d^K^{E), includes the ef- 
fects of nonlocal self-energies. Equation ([9]) is exact in 
the atomic limit (Appendix |A|) . and is a good approx- 
imation for large disorder, where the imaginary part of 
5e {E) is small and independent of E (see Appendix |B|) . 
equation is derived assuming U < Uc, where C/c ~ A 
in the large disorder limit, since the system is a gapped 
Mott insulator for U > Uc- 

Equation ([9|) gives an explicit relation between p{E) 
and the functions Se(a;) and Ae(w). One can interpret 
the derivatives 9eEe(a;) and 3eAe(a;) as the response of 
the self-energy and hybridization function to changes in 
the local potential or, equivalently, the response of these 
functions to the disorder potential. This is reminiscent of 
the situation in weakly correlated metals, where a similar 
analysis related the ZBA to the response of the charge 
density to the disorder potential 

We use Eq. ([9]), in conjunction with our numerical cal- 
culations, to establish the relative importance of Ai(cij) 
and T,ii{uj) in forming the ZBA. The first step is to ex- 
tract 9eEe(aj) and dcA^{Lu) from numerics. This process 
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is illustrated in Fig. [2l For a given disorder configura- 
tion, Y,ii{E) is calculated from Eq. at a fixed value 
of E, chosen to be ep in Fig. [2j (Note that, for a fi- 
nite size lattice, both Ai(uj) and Sij(aj) are real.) The 
collected values of 'Sii{E) for all sites i and for 1000 con- 
figurations are shown in Fig. [2l[a). Data is shown for 
the two ranges, ei + U ~ sp and ~ sp^ that con- 
tribute to piep)- A disorder-averaged ^e{E) is found 
by making least-squares quadratic fits to the data in 
each range, from which the derivatives d^T,^{ep)^=^f, and 
d^T,^{ep)^^sj,-u are extracted. An identical set of calcu- 
lations is then made for deA^{£p). The calculations are 
repeated for other values of E, and resulting derivatives 
are plotted as functions of E in Fig.[2][b). 

As a check, we compare in Fig. ^c) the DOS from 
Eq. ^ , calculated using the values shown in Fig. [2][b) , 
with the exact DOS. The agreement between the two 
is very good. We have repeated this analysis for other 
values of A, and continue to find qualitative agreement 
down to the Mott transition at A « [/ [Fig. [2];d)]. 

Figure ^h) shows the relative contributions to the 
ZBA made by T,ii{uj) and Ai{ijj). The figure shows that 
dj:^{E)^^-^ is negative for both LHO [E = E) and UHO 
{E = E—U). From Eq. ([9]), we see that a negative deriva- 
tive corresponds to an increase in p{E), and not to the V- 
shaped suppression of the DOS required to form a ZBA. 
This result demonstrates that the ZBA does not come 
from the local self-energy, and is therefore not a remnant 
of the Mott gap. More significantly, it demonstrates that 
the physics underlying the ZBA cannot be reproduced 
by approximations that include only the local self-energy, 
such as single-site DMFT or the Hartree-Fock approxi- 
mation. The ZBA that appears in unrestricted Hartree- 
Fock calculations must have a different origin than that 
found here. 

In contrast to the self-energy derivative, d^A^{E) ^^-^ is 
positive for E Ki ep and negative away from ep^ indicat- 
ing that interorbital hybridization shifts spectral weight 
away from ep. This shows that the ZBA comes from non- 
local correlations embedded in the hybridization func- 
tion. On the one hand, this is not surprising since the 
Hubbard model in low dimensions is known to map onto 
effective models with nonlocal interactions; on the other 
hand, the energy scale t of the ZBA is not consistent with 
the energy scale /U of these effective models. 

We note one further interesting feature of Fig. [Hb): 
the plots of deA^{E)f=E and deAt{E)f=E-u are asym- 
metric with respect to ep. This asymmetry indicates 
that LHOs and UHOs behave differently when they are 
below or above ep. We will return to this point below. 

In summary, we have established two main results in 
this section. First, we have developed an expression, 
Eq. dni), for the DOS that relates p{E) to the response of 
^e{E) and Af{E) to the disorder potential. Second, we 
have used this expression to analyze exact diagonaliza- 
tion results, and have shown that the ZBA is the result 
of nonlocal correlations, rather than the local self-energy. 



B. Structure of the Hybridization Function 

In the previous section, we established that the ZBA 
can be related to the derivative of A^{ui) with respect to 
the site energy e. In this section, we analyse the structure 
of Af{uj) in more detail in order to see the role of spin 
and charge fluctuations in forming the ZBA. 

We begin by writing Ai(uj) in terms of an alternative 
exact expression's that is more transparent than the orig- 
inal definition [Eq. (|4])]: 

A,{UJ) = + + ^M], (10) 

where G^j,(aj) is a Green's function matrix element for 
the lattice with site i removed4S This equation shows 
explicitly how the matrix elements Uj + Sy (aj) couple 

the site i to the rest of the lattice. In general, G^j.(cj) 
is not trivial to calculate, and this expression is of use 
only when G^^(a;) can be simplified through some ap- 
proximation or limit. Here, we are in the limit of large 
disorder and low dimension, for which G^-^{uj) is approx- 
imately local. In our discussion, we thus consider only 
the dominant contributions, with j = fc, in the sum in 
Eq. ([101); 

A,{lo)^ ^ [-t + i:,,{u)fG%{uj), (11) 

where j G nn^ indicates that j is a nearest neighbor of i. 

We note that, while Ai(w) is real for a single disor- 
der configuration on a finite lattice, the disorder-averaged 
hybridization function Af{ui) is complex. The real part 
of Afiuj) describes shifts of the LHO and UHO ener- 
gies while the imaginary part describes the broadening of 
these orbitals due to the lattice. For the analysis in this 
work to make sense, the broadening must be much less 
than the level spacing (~ A/22;) of the local spectrum, 
so that discrete energy levels at each site keep their dis- 
tinct identity. We can estimate the broadening from a 
simplified disorder average of Eq. ([TT|). Setting = 0, 
we obtain 

A° (a.) = zt^{G%{Lo)),, (12) 

where the sum over j is replaced by the factor z, and 
(. . .)j = A~^ ■ • ■ ^^i the disorder average over 

site j. This equation assumes that G^jj{uj) with different 
j are independent of each other. The imaginary part of 
Eq. (USD is 

ImA°(a.) « -.ziV-)--^, (13) 

which gives a broadening of 0(zi^/A). The condition 
that this is much less than the level spacing of the local 
spectrum can be written 2z^t^/A^ <C 1, which is met 
provided our initial assumption 2zt/ A ^ 1 is met. 
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Equation shows that the nonlocal self-energy is 
central to the ZBA. To proceed further, we need an an- 
alytic form for this self-energy, and we adopt a partial 
fractions expansion for the self-energy that is based on 
the equation-of-niotion method.'^ The rationale for this 
choice is that the equation-of-motion method correctly 
reproduces the LHO and UHO in the atomic limit, and 
has been shown to be accurate for the two-site AHMM' 
In general, we expect this method to work well when 
short-range physics dominates. The nonlocal self energy 
has the form 



(W - eicr - Uhiff){uj - eja - Uhja) ~ - ... 

"(14) 

where we suppress the explicit dependence of (w) and 
Pij on cr because we are considering only nonmagnetic 
phases, where hi-^ = 1 — 71^^ — {fiia), with a = —a, 
and where 

p,, = (Sn^^Shj^) + {S,+ S,.) - (dId,). (15) 

(Here, (. . .) indicates the expectation value, rather than 
the disorder average.) The three nonlocal correla- 
tions making up pij involve density fluctuation operators 
Sfii^ = fiia- — riia-, spin- flip operators Si±, and pair anni- 
hilation operators Di = CiiCi-\-. The last of these three is 
an order of magnitude smaller than the other terms and 
is discarded for the remaining discussion. 

In general, the usefulness of Eq. (|14p is limited by the 
difficulty of finding the higher-order terms in the contin- 
ued fraction. These terms are important for determining 
the pole structure of the self-energy, but do not change 
the fact that E^- cx pij. In the disorder- free Hubbard 
model, it has been shown that these higher order terms 
are qualitatively important*^ however, the strongly dis- 
ordered case is close to the atomic limit and may be un- 
derstood qualitatively through a truncated self-energy, 
obtained by dropping the 0{t^) term in (fT4|). We check 
this assertion numerically: we calculate an approximate 
Ai{uj) using the self-energy pi)) in Eq. pT|). and then cal- 
culate an approximate DOS using Eq. ([9|). The results 
are plotted in Fig. ^c) in comparison with exact diag- 
onalization calculations, and the agreement between the 
two is good. 

We showed in the previous section that the ZBA comes 
from the response of (E) to the disorder potential via 
the derivative deAe{E). The main idea suggested by 
Eqs. pT|) and (|14l) is that this response is directly re- 
lated to the response of j (E) , and therefore of pij , to 
the disorder potential. We show in the next section that 
there are other contributions, but that a large part of 
the ZBA can indeed be traced back to the response of 
the nonlocal charge and spin correlation functions to the 
disorder potential. 

We note that the form of Ey (w) explains the asym- 
metry in deA^{E)f=E and d^Af (E)^=E~u with respect to 



£F^ shown in Fig. [2Kb). This figure shows that the ZBA is 
formed from a shift away from ep oi LHOs below e p and 
of UHOs above ep- According to Eq. (ITU) , this asym- 
metric shift occurs because the correlation pij is largest 
when sites i and j are both singly-occupied, namely when 
ep — U ^ £ij £j ^ £f- (The spin correlations vanish when 
either site is empty or doubly occupied.) This condition 
on €i and tj is equivalent to the requirement, at each site, 
that the LHO be below ep and the UHO be above ep. 

In summary, we have used a form for the hybridization 
function that shows explicitly the role of the nonlocal 
self-energy. We have proposed using an analytic form, 
Eq. for this self-energy, and have shown numerically 
that it reproduces the density of states obtained by exact 
diagonalization. The main result of this section is that 
the nonlocal self-energy, and therefore the ZBA, depends 
on nonlocal spin and charge correlations. 

Ideally, one would now like to use this formalism to de- 
rive an analytic expression for the density of states; this 
requires knowledge of pij and is in general quite difficult 
since pij is different along every bond in the lattice. In 
the next section, we therefore focus on a simple model 
for which pij is known, and the DOS can be found ana- 
lytically. 



C. Density of States 

As a simple application of the formalism derived in the 
previous sections, we calculate the DOS for the two-site 
AHM (2SAHM). This model has been studied elsewhere 
by direct diagonalization of the Hamiltonian,-^"'^^ and 
provides a point of comparison for the current work. Our 
approach is straightforward: we use the self-energy pi)) 
to find an approximate hybridization function with which 
we evaluate the density of states using Eq. ([9]) . 

The 2SAHM consists of an ensemble of two-atom 
"molecules" with random site energies and tj. The 
disorder averaged hybridization function for site i is 

In this form, the hybridization function has a useful sym- 
metry f Appendix (P)) 

5AuHo(£F + E) = aALHo(£F - E), (17) 

where 9Alho(-E') and 9Auho(-£') have similar definitions 
as 9S'lho(£') and dS\jYio{E), and E is the energy E mea- 
sured relative to e^. 

One consequence of this symmetry is that contribu- 
tions to 9Alho(-E') that are even under E ~E are 
more important for the ZBA than those which are odd. 
To show this, we define Sp{E) to be the change in the 
DOS due to the hybridization function, namely 6p{E) = 
p{E) — po{E), where po{E) is evaluated with d^A^^E) set 
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to zero. To linear order in d^A^^E), Eq. ([9]) gives 



5p{E) 



1 



[l + 9I]LHo(i?F [l + 9SuHo(i?P 

(18) 



Noting, from Fig. d^b), that 9Slho(S) ~ d^^BoiE) 
near ei?, we get 



SpiE) 



1 aALHo(^) + 9AuHo(^) 
A [l + 9S](ei.)]2 



(19) 



From this, and from Eq. p7)) . it fohows that the most 
significant contributions to 5p{E) come from terms in 
9Alho(-£') that are even in E. 

To calculate 9Alho(-E'), we expand Eq. as 
A,(cj) = AO(w) + A^(w) + A^'(lj), where 



K^{u:) - -2t(S,,(a.)G5^.(a.)), 



= (E.,(a.)2G;^.(c.)), 



(20) 
(21) 
(22) 



We have evaluated each of these terms analytically and 
find that, by far, the largest contribution to the ZBA 
comes from 9Alho(-^)- particular, 9Alho(-^) 
odd function of E and therefore makes almost no con- 
tribution to the ZBA; 9Alho(-^) contains both odd and 
even terms and therefore does contribute to the ZBA, 
but is an order of magnitude smaller than dA'^^Q{E). It 
is, perhaps, not surprising that the term containing the 
highest power of (oj) makes the largest contribution to 
the ZBA. For clarity, we include only results for A'l{E) 
in our calculation of 6p(E). 



Using Eq. for G^.Aio), we obtain 



A/2 



^ J~A/2 ' {UJ ~ ^t~Uhi-jY{iO 

and differentiating this with respect to e^, we obtain 



uy 



(23) 



To calculate OA'^^q {E) , we set i 



A/2 



1 



- Uh 



Uhj-^){uj - ej)(w ■ 



U) 



(24) 



E in Eq. ^IQ. Then there are four terms, proportional to d^-pf^^ to pf^de^hi-ff, 



to p^^de^hj-ff, and to -. The last of these is a factor t/U smaller than the others and is discarded. 



Because of the simplicity of the 2SAHM, we can write 
the coefficients pij, hi^, and hj-^ in terms of the many- 
body wavefunction for the two site system, and thus find 
their explicit dependence on and ej. This makes the 
integration over ej possible. The calculations are com- 
plicated by the fact that we do the integration at fixed 
chemical potential, meaning that the number of electrons 
in the ground state depends on Ci and ej. The domi- 
nant contribution to the ZBA comes from cases where 
the ground state has two electrons, and we include only 
this term in our result. The calculations are somewhat 
lengthy, and we leave the details to Appendix [Cl 

The result of these calculations is, from Eq. ([T5| and 
Eq. dClll), 



6p{E) 



-28V2t 



27A2(l + aE)2 



F2{x)+F,{x) + - 



(25) 



where x = {2t'^ ~ E'^)/ {2V2t\E\), with E ^ E ~ sp and 



^2(2;) 



+ tan ^{x), 



F,{x) = -F,{x)-j^^^^. 

Equation ([25]) is plotted in Fig. |3] for the case A — 20i. 
For this plot, the unknown prefactor 1 -I- 91] is taken to 
be 0.7, based on the value of 9eEe(£')|e=e^ in Fig. [S^b). 
The resulting plot is qualitatively consistent with exact 
resuhs for the 2SAHM35-36; from Eq. the width of 

the ZBA is of order 2\/2t^ and the depth is proportional 
to t//\^. 

In previous studies of the 2SAHM, the ZBA was at- 
tributed to level repulsion between many-body states. 
Here, level repulsion is implicit in 9Alho(£') f^nd 
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FIG. 3: Theoretical density of states from Eq. (|25|) . Results 
are shown for A/t = 20. Note that the curve is independent 
of U. 



<9AuHo(-E')i since these describe the shifts of the atomic 
LHO and UHO due to neighboring sites. These shifts 
are primarily due to the response d^. (E) of the non- 
local self-energy to the disorder potential. In Eq. (1^ . 
we showed that di-Tjij{E) depends on the local charge 
susceptibilities de^hi^ and d^hj-ff, and a generalized sus- 
ceptibility deiPij. For the 2SAIIM, the last term is the 
largest, so that the ZBA is mostly due to the response of 
the nonlocal spin and charge correlation functions mak- 
ing up Pij. 

It is interesting to note that Mott physics suppresses 
this response. This is because the local Coulomb interac- 
tion tends to fix the charge density at each site, so that 
the spin and charge correlations are only weak functions 
of ei. For example, configurations in which and ej are 
near ep have a singlet ground state |s), with corrections 
of order t/U. A small change in changes this ground 
state, and therefore pij, by order t/U. Thus d^-pij is 
suppressed by Mott physics. This is not the case when 
Ci and €j + U are within t of ep- Then \s) and |02) are 
nearly degenerate, and the proportions of |s) and |02) 
making up the ground state vary linearly with . In this 
regime, d^Pij is not small. The ZBA therefore comes 
from disorder configurations in which Mott physics does 
not suppress nonlocal charge fluctuations. 

The results presented in this section are valid for 
A > t/ ^ t. When C/ > A, the spectrum has distinct 
lower and upper Hubbard bands. In our calculations for 
the 2SAIIM, the ZBA collapses rapidly when the Hub- 
bard bands no longer overlap, since conflgurations with 
degenerate LHO and UHO no longer occur. This appears 
to contradict results reported by Chiesa et al.,^^ where 
the ZBA persisted for U > A, away from half-filling. Di- 
rect comparison with Ref. [l^ is not straightforward since 
they are not in the regime A ^ 2:t in which our theory 
is valid. We have performed preliminary exact diagonal- 



ization calculations for one- and two-dimensional clusters 
for the case U > A ^ zt; these show that while the slope 
of the ZBA (namely, dEp{E)) is approximately indepen- 
dent of [/, the width and depth are stronger functions of 
A than when U < A. We find that the width of the 
ZBA is not simply t in the gapped phase; however, these 
results are preliminary, and a careful study is required to 
resolve this discrepancy. 

III. CONCLUSIONS 

In this work, we have discussed the origins of the 
disorder-induced zero bias anomaly in the Anderson- 
Hubbard Model. Several aspects of this zero bias 
anomaly are unique to strongly correlated systems with 
short range interactions. Most significant is the fact that 
the width of the anomaly is set by the hopping matrix 
element i, and is independent of the interaction strength 
U and disorder potential A over a wide range of A and 
U . In the two-site Anderson-Hubbard model, this has 
been understood as the result of level repulsion between 
lower and upper Hubbard orbitals<^i2^ 

Here, we have gone beyond the 2SAHM, and have 
shown that the underlying physics of the zero bias 
anomaly in larger clusters can be extracted from an anal- 
ysis of exact diagonalization calculations. The analy- 
sis is based on an expansion around the atomic limit, 
and is appropriate for disorder A much larger than the 
clean-limit bandwidth zt. Through this analysis, we have 
found that the local Coulomb interaction generates non- 
local spin and charge correlations between adjacent lat- 
tice sites, which cause an overall shift of spectral weight 
away from the Fermi energy ep- By this mechanism, a 
V-shaped zero bias anomaly is formed in the density of 
states aX £p. 

Specifically, the zero bias anomaly comes primarily 
from the response de/Tjij{E) of the nonlocal self-energy 
to the disorder potential. Mott physics tends to suppress 
this response; however, disorder configurations in which 
many-body Fock states are nearly degenerate are sen- 
sitive to small changes in the lattice potential, and for 
these configurations dcJ^ij{E) is not small. 

Using the formalism developed in this work, we have 
obtained an analytic expression for the DOS of a two-site 
Anderson-Hubbard model. This expression reproduces 
the essential physics of the zero bias anomaly found nu- 
merically; the anomaly has a width of order t, and a 
depth which is independent of U when U ^ t. 
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Appendix A: Results for the Atomic Limit 

In the atomic limit, the exact local Green's function in 
the nonmagnetic phase is 

1 



cj - - Sii(w) 
where the exact self-energy is,— 

E,,(^) -U- + 

and the charge density is 



2, ei + U <eF 

1, £f — U < ei < Ef 
ei> Ef 



(Al) 



(A2) 



First, we note that it follows directly from (|Aip and (|A2p 
that the term S^{E) — U in Eq. ([8]) vanishes identically 
in the atomic limit. 

Next, we check that Eq. (|9]) is exact in the atomic limit. 
From Eq. 



n,/2 



(A3) 
(A4) 



'^•=^-^ - n,/2 
Taking, for example, E slightly less than e^? we obtain 



piE) 



1 

A 
1 

A 
3 

2A' 



1 



i + a,E,(E) 
1 1 



+ 



1 



e=E-U 



1+1 1+0 



(A5) 



for the disorder-averaged DOS. This is the exact result. 



where the argument of the logarithm is complex. E has 
an infinitessimal positive imaginary part so that Eq. (jB2p 
reduces to Eq. © when S and dS are real. 

In this work, S^{E) is complex as a result of the dis- 
order averaging process. We find that the hybridization 
function introduces imaginary components. 



S S-iT 
dS dS + i-f 



(B3) 
(B4) 



where F ~ zt^/A and 7 ~ zt/A. Near the Fermi energy 
at half filling, E ~ U/2 such that 



(±^-£;)(i + 9^) 



(4-.) 



> \s-ui 
1 » r, 



except near the Mott transition at [/ « A. Then 



P{E) 



- y 

ttA _ ^ 1 1 + as* 

B=B.B-f/ 



tan 



E=E,E-U 
1 -7 



tan 



-1 7 



l + dS 



-{l + dS) 



7 



- y 



l + dS 

7 



In 



i^-E) 



(t + e) 



ttA _ ^ \l + dS l + dS 

E=E,E-U 



In 



E 



E 



(B5) 



where tan-i [-7/- (1 + 95')] w -7T + j/{l + dS). The first 
term in Eq. (|B5I) is the result found in Eq. ([9|) , while the 
second term increases p{E) by order ztj A? . This term is 
comparable in magnitude to the corrections responsible 
for the ZBA, but is featureless near E = ef, and therefore 
does not contribute to the ZBA. The conclusion to be 
drawn from this appendix is that the expression Q is 
sufficient to understand the ZBA provided zt/A<^ 1. 



Appendix B: Density of States for Complex Self 
Energies 

If the self-energy S^{E) is complex, then the analysis 
leading to Eq. ^ is unchanged, 



G,iE) 



1 



{e-E)[l + dS] + [S-U]' 



(Bl) 



where we use the compact notation dS — diS^{E)\^^-^ 
and S = S-j^^E); however, the disorder-averaged density 
of states is 



PiE) 



E=E,E~U 



In 



{^^E){l + dS) + S-U 
{-%~E){l + dS) + S-U 



(B2) 



Appendix C: Derivation of 9Alho('^) 

Here, we calculate the derivative dK'l^^o (E) for an en- 
semble of pairs (i, j) of isolated sites with random ener- 
gies. The Green's function for site j with site i removed 
is the atomic Green's function 



bJ — €j UJ — ej — U 

uj — €j — Uhj-ff 
{uj - ej){uj - - Uy 



(CI) 
(C2) 



where cr = — cr, /ij^ = 1 — rij-^, and where we suppress 
the spin index on G in the nonmangnetic state. 
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FIG. 4: Phase diagram for an isolated pair of sites 

with site energies ei and ej. The figure shows the number 
Ne of electrons in the ground state, and is divided into four 
quadrants by the dashed blue lines. The quadrants are la- 
belled A,. . . ,D. For the 2-electron ground state, the most 
important regions are ei ^ ej + [/ ~ £f (quadrant D) and 
ti + U ~ Ej ~ £f (quadrant A), which correspond to the LHO 
on one site and the UHO on the other being nearly degenerate 
with ef- For the 1-electron and 3-electron ground states, the 
most important regions of the phase diagram are ei ~ ej ~ £f 
(quadrant B) and ei ~ e^ ^ ef — U (quadrant C) respectively. 



Using Eq. del 



and 



KM 



eu^ /•A/2 



de 



A/2 (w - Ei - Uhi-sY{LO - ej - Uhj-^){LO - ej){i^ - tj - U) ' 



(C3) 



5A'/.(w) 



2 _ 
-1 



Uht. 



+ Uhj^ - ei){ej - ei){ej + U - et) 



(C4) 



As discussed in the main text, the term in the square brackets proportional to pfj/Uhiff is a factor t/U smaller than 
the other terms, and is discarded. 



Each pair of sites in the ensemble may have anywhere 
from to 4 electrons, depending on and ej, and in 
order to evaluate the integral in Eq. (jC4p . we need to 
keep track of the different states. Figure |3] shows that 
there are four different possible ground states when 
is fixed near ep, having a total of 0,1, 2, or 3 electrons 
shared between i and j. 

The 0-electron ground state does not contribute to 
9Alho(^') because pij — 0. The 3-electron case also 
does not make a substantial contribution, in this case 
because the derivatives d^.hf^^ d^-hjs^, and de^Pij are of 



order t^/[/. This follows because, in region D of the phase 
diagram (Fig. the ground state wavefunction has the 
form 

|3e)«|(72) —\2a), 

e.. - ej 

where |(t2) indicates that site i has a single spin-a elec- 
tron and that site j is doubly occupied. Because — Sj ^ 
U, the derivatives in Eq. (jC4[) are of order t/U'^. 

The remaining contributions to Eq. (jC4[) are from the 
1-electron and 2-electron ground states. It turns out that 
is an order of magnitude smaller in the 1-electron case. 
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where (Si+Sj-) = 0, than in the 2-electron case, and so 
we focus our attention on the latter. 

For ei ^ Ep, the most important contributions to the 
ZBA for the 2-electron case come from quadrant D of the 
phase diagram in Fig. |4l In this quadrant, there are two 
important configurations: the singlet \s) = (| — | |t 
))/\/2, and the double-occupancy state |02) which has 
both electrons on site j. These states have nearly the 
same energy, so we use degenerate perturbation theory. 
We project the AHM onto \s) and |02) to get the Hamil- 
tonian matrix. 



2e 



-V2t 



(C5) 



-%/2t 2ej+U 

from which follows the ground state wavefunction, |2e) — 
a\s) + VI - a2|02), with 



a 



(C6) 



where y = {ej + U — ei)/2. We can write the expectation 
values in Eq. (|C4I) in terms of a^: 



a 
Y 



P^3 



'JO- 

1-'- 

4 



a 

T 



(C7) 
(C8) 



It simplifies our calculations significantly that the deriva- 
tives in Eq. ()C4p all reduce to derivatives of . We use 



■^dy and substitute de, — > 2dy to obtain 



dA'l (uj) 



yo 



(C9) 



We need the principal part of the integral for SA^jjq (E) . 
In deriving this expression, we have neglected terms of 
order t/U. 

The integration limits are given by the range of y over 
which the ground state has two electrons. It can be 
shown that the 2-electron state is stable in region D of 
Fig.Hfor^ 



(CIO) 



where ii 



Ep- The boundaries of the 2-electron 



phase in this region are shown as thick black lines in the 
figure. Setting ei = E, we obtain the integration limits 



E 



yo 



E 2 



yi = oo, 



E 

"^'^E-r 



[E < 0) (Cll) 
(^>0), (C12) 



where E — E — Ep. The integration limits at ±oo come 
from the boundaries of region D, which are taken to be 
far from Ep and Ep — U. This assumption does not change 



our results significantly because the integrand in Eq. (IC9[) 
is peaked near y = 0, because of the factor 

d - I. 

" ~ (y2+2t2)3/2- 

We, for the same reason, can expand 

-4 / . \ 4 / 



1 



1 



y 



3^y^ + 2t\ 
4y 



10y2 



3^y^ + 2^2 9(y2 + 2i2) 



and 



1 



a 
T 



to obtain 
dA'l^^oiE) - Re 



8 8VF+2t2 



dei 



4^/2t 
9A 



7 



F, 



"1 VI 

yo 
(C13) 



where F^{x) = n / x"" V(l + a;^)^^"^^- The func- 
tions Fn{x) are odd (even) when n is even (odd). Be- 
cause the limits yo and j/i are odd in E, F2{x) and -F4(x) 
actually make a contribution to (9Alho(-E') that is even 
in E. Using the symmetry of Eq. P7|) . it follows that 



28V2t 



\V2tJ 



27A 



yo 



Explicitly, 



F,{x) 
Fi{x) 



x^ 
1 



+ tan ^{x). 



(C14) 

(C15) 
(C16) 



Appendix D: Symmetries of dAi^uo{E) and dAuuo{E) 

In this appendix, we prove the relation 9Auho(-E') — 
9Alho (—-£') (for convenience, we take Ep = in this 
section). This result is based on the symmetries of the 
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single-site Green's function, Eq. (jCll) , and the self-energy 
Eq. ((HD. 

The proof proceeds as follows: 9Alho(-E') has contri- 
butions from 2-electron states in region D of Fig. 0] and 
1-electron states in region B; 9Auho(^') has contribu- 
tions from 2-electron states in region A of Fig. |4] and 
3-electron states in region C. We show that there is a 
correspondence between regions A and D, and between 
regions B and C, with the result that d^K^{E)\^=E in B 
(or D) is equal to d^h-^{—E)\^^E-u in C (or A). 

Suppressing subscripts, we write d^K^{E), with h-e{E) 
given by Eq. p^ . as 

dk = 2{{-t + S)(ai])G) + (i-t + T.fdG). (Dl) 

Now consider a pair of sites with and tj belonging to 
region D, and a corresponding pair of sites with = 
— U, €j — ej + U belonging to region A. For region 

D, the wavefunction is |2e) — ay\s) + ^1 — a^|02) with 
y — {ej + U — ei)/2, while for region A, |2e)' = ay\s) + 
^1 — a^|20) with y = {ei + U — ej)/2, where ay is the 
same in both cases. Because of this symmetry, h'^-^ = nj-^ 
and n'j-^ = hj-ff. It follows immediately that the local 
Green's function Eq. (jCl|) satisfies 

G^, = -Ge^=E, dG'^.^u^E = dGei=E (D2) 

where G' is the Green's function for primed site energies, 
and G is for unprimed site energies. It also follows that 



Sij(i?) — >■ T,{y) for regions A and D, with S(y) the same 
even function of y in both cases, but with y specific to 
each region, as above. Thus 

^'ei=E-U = ^ei=E, dY,'ei=E~U — —9^ei=E- (D3) 

Equations (|D2I) and (jD3p suggest that dA is even under 
(ei,ej) — ^ (e^,ej); however, an additional negative sign 
arises from averaging over j. For e^. 




while for e'^ 




Because of the inverted integration limits, we obtain 
(considering only contributions from regions A and D), 

d,A,{E)U=E = d,A,{-E)U^E-u- (D4) 

An identical result is found if we consider primed and 
unprimed site energies belonging to regions B and C re- 
spectively, which proves 

dAi^noiE) = dAu^o{-E). (D5) 
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